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AN APL FUNCTION FOR BIVARIATE NORMAL PROBABILITIES 



INTRODUCTION: 


Bivariate 


normal 


distributions have 


many 


applications 


such as. in 


combat modelling. 


weapons 


systems 


effectiveness 


studies and 


weather prediction 


problems ; 


several 


other applications are discussed in 


[4]. It 


is well known that 


probability 


statements 


for a 


general 


bivariate 


normal 


distribution 


can be transformed into 


equivalent statements 


for a 



Standard bivariate normal distribution with 0 means and variances 
equal to 1. It is thus sufficient to be able to compute 

cumulative probabilities for a standard bivariate normal 
distribution with probability density function 

f(x,y) = „ L 2" e _oo<x<°°-°°<y<°° 

2 TTvl - P 

Similar to the case with the univariate normal distribution, the 
bivariate cumulative distribution function (c.d.f.) 



F(h,)c) 



P[X<h,Y<k] 



k h 

f(x,y) dx dy 

—00 —00 




does not have a closed form solution and numerical integration is 
the usual approach to evaluating such integrals. Because of the 
importance of the distribution, the National Bureau of Standards 
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(NBS) has published an extensive set of tables for P[X>h,Y>k] 

[5] . Another set of tables using a different approximation, has 
been generated by Owen [6]. Other ways to approximate the 
bivariate normal integral are discussed in [1], [3] and [7]. All 

of these approximations involve numerical integration and require 
extensive computer programming, rendering them to be not readily 
suitable for obtaining on the spot results. These days, with the 
easy availability of microcomputers, it would be useful to have a 
program to compute bivariate normal probabilities interactively 
and also be able to incorporate such a program within other 
application programs. This will allow an analyst to perfomm 
sensitivity studies by varying the input parameters in an 
application program and observing the effect on the measures of 
effectiveness of interest. 

Recently, .Wang [8] developed a new algorithm to compute 
bivariate normal probabilities, that does not require numerical 
integration, relatively easy to program and provides quite 
accurate results. Investigations by Wang indicate that the 
computed probabilities compare very well with those in the NBS 
tables and the computer resources needed are not excessive. The 
approach used by Wang is to start with an approximate contingency 
table of cell probabilities for a rectangular grid on the x, y 
plane and then apply the '‘iterative proportional fitting 
algorithm" [3; sec 3.5] to modify the cell probabilities; the 
iteration is continued until the marginal probabilities in the 
contingency table coincide with their true values to within a 
specified degree of accuracy. Since the marginal distributions 
of the random variables X and Y are univariate normal, the exact 
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marginal probabilities can be determined using computer programs 
available in most statistical software packages or by writing a 
subprogram for the purpose. It should be noted that, although a 
bivariate normal distribution is defined over the entire plane, 
for all practical purposes it is sufficient to consider only the 
square subregion [-4,4] x [-4,4] since the probability content 
outside of this subregion is negligible. 

Wang's algorithm is defined by the following steps: 

1. Let -4 = aQ a^^ ... ajj^ = 4 be a partition of 

the interval [-4,4] along the x - axis and -4 = bQ b^^ 

... bj^ = 4, a partition along the y - axis; identical 
equal length partitions for both axes is recommended to 
reduce computational complexity. 



2. Let 


Ax = Xj^ - 3 - 2.-1 


i = 1, 2, . . . , 


m 




Ay = bj - bj_i 


j = 1, 2,..., 


n 




iHl = + a^ 


i = 1, 2 . . . ,m 






nj = bj_i + bj 


j = 1,2, ... ,n 




where 


p = p (1 - Ax^) (1 - Ay^ ) 

12 12 

is the correlation coefficient 


(specified) 


3 . Let 


= P ^ 1 


X L^^)***) in 





p. = P [b. ,< Y < b.] j = 1,2,.. .,n 

J J-1 - J 

be the marginal probability contents of the i-th and j-th 
subintervals along the x and y axes respectively. 
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4. Compute pj^j , the starting approximate probability of the 
(i, j) th cell as 




P m.n. 
1 J 




i = 1,2, ...m; j = 1,2, ...,n 



5. The application of the iterative proportional fitting 
algorithm results in the following equations for the modified 
probability of the (i,j) th cell after the k-th iteration: 



6 . 



,(k) 

ij 



(k-1) 
,(k-l) 

ij 



E 



j=i 



k = 1,3,5,... 



k = 2,4,6, .. . 



(k-1) 

P. . P. 
iJ J 

I 1=1 

Continue the iteration process until for some 



n 



E p ' - 5, 

j=l 



< £ 



and 



m 



E p 

i=l 



(k) _ p 

ij j 



even number k 

z 



where e is a prespecified degree of accuracy with which the true 
marginal probabilities agree with the marginals in the 
contingency table. 



7 . To compute P[X£h, Y £k] (orP[X>h, Y> k] ) sum 

the probabilities in the contingency table over those cells for 

which £ h (a^ > h) and < k (b^ > k) . In those cases where either 

h 4 a and/or k b . for any of the partition points a . and b . , the 
i J 
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accuracy of the approximation can be improved by including h 
and/or k as additional partition points. 

An APL function, called BVN, to compute bivariate normal 
probabilities (both P[X<h, Y<k] and P[X>h, Y>k]) is 
presented in the appendix. This function invokes another APL 
function called NCDF to compute the marginal univariate normal 
probabilities. The BVN function can be run on an IBM- PC / AT 
compatible microcomputer using an APL language system such as 
APL* PLUS from the STSC corporation; the function can also be run 
under VSAPL on the NPS mainframe computer. The function runs 
interactively and calls for keyboard input of the desired equal 
length partition size for the x and y axes and the degree of 
accuracy e in approximating the marginal cell probabilities. 
With only minor modifications, the function can be imbedded 
within another APL function as a subprogram. Computations using 

the BVN function indicate that with partition size x = y = .2 

for both the x and y axes and e = .00005, a four decimal place 

accuracy as compared with the NBS tables can be achieved. The 

computational time, as is .to be expected, increases with a 
decrease in the partition size x or y, an increase in and 

to a lesser degree a decrease in e . With x = y = . 2 , 

e = .00005 and .1< <.8 the computational time (clock time) 

was between 90 and 180 seconds on the Zenith Z-248 (an AT type) 
microcomputer, and on the IBM 3033 mainframe computer these times 
were between 1 and 20 seconds. The computational time can be 
reduced considerably (to about 30 seconds on the Z-248) by 
choosing x = y = .5 but then only a three decimal computational 
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accuracy can be expected. For a fixed value of if the c.d.f. 
is to be calculated for several choices of (h,k), the BVN 
function needs to be run only once; with a very minor 
modification the contingency table of bivariate cell 
probabilities can be saved in a matrix and all that is left is to 
sum the probabilities of the appropriate cells. If needed, it 
would be quite straight forward to generate tables for various 
choices of and (h,k) . 

Professor I. O 'muircheartaigh of the O.R. department and I 
are in the process of completing a Fortran program for the 
problem and expect to submit the code for publication. 
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APPENDIX 



THIS APPENDIX CONTAINS THE LISTING OF TWO APL VARIABLES BVNHOW AND 
NCDFHOW AND TWO APL FUNCTIONS BVN AND NCDF. THE TWO HOW VARIABLES 
PROVIDE SHORT DESCRIPTIONS OF THE COMPUTATIONAL SCHEMES, THE INPUT 
PARAMETERS AND THE SYNTAX FOR THE TWO FUNCTIONS. 

BVNHOW 



THE FUNCTION BVN COMPUTES THE C.D.F. OF A STANDARD BIVARIATE 
MORMAL DISTRIBUTION WITH CORRELATION COEFFICIENT USING AN 
(ALGORITHM PROPOSED BY YUCHUNG J. WANG (BIOMETRIKA, 1987, NO. 74, 185-90). 
THE ALGORITHM CONSISTS OF PARTITIONING THE X-Y PLANE INTO RECTANGULAR 

:ells, an initial approximation of the cell probabilities and an 

ITERATIVE SCHEME TO MODIFY THE CELL PROBABILITIES. THE ITERATION 
PROCESS IS TERMINATED AS SOON AS THE MARGINAL PROBABILITIES COINCIDE 
«JITH THEIR EXACT VALUES (THAT ARE UNIVARIATE STANDARD NORMAL 
>ROBAB I LI TIES, COMPUTABLE USING THE APL FUNCTION NCDF) TO WITHIN A 
5PECIFIED DEGREE OF ACCURACY e . THE SYNTAX FOR THE FUNCTION IS 

I RHO BVN W 

IHERE RHO IS THE CORRELATION COEFFICIENT AND W = <x,y) IS THE POINT 
)T WHICH THE C.D.F. IS TO BE COMPUTED. THE FUNCTION WILL CALL FOR 
i'HE INPUT OF THE DESIRED LENGTH FOR AN EQUAL PARTITIONING OF THE 
NTERVAL C"4,43 ALONG THE X AND Y AXES AND e THE DESIRED DEGREE 
F ACCURACY IN THE MARGINAL PROBABILITIES. THE RECOMMENDED CHOICES 
RE .2 FOR THE PARTITION LENGTH , AND .00005 FOR e . THE TOTAL 
□MPUTATION TIME, ON THE ZENITH Z-248 MICROCOMPUTER, SHOULD BE BETWEEN 
0 AND 180 SECONDS DEPENDING ON THE SIZE OF RHO. 

NCDFHOW 



HE FUNCTION NCDF COMPUTES THE C.D.F. OF A STANDARD NORMAL DISTRIBUTION, 
SING THE APPROXIMATION DEFINED IN EQUATION 26.2.17, PAGE 932 IN 
HE HANDBOOK OF MATHEMATICAL FUNCTIONS, M.ABRAMOWITZ AND I.A.STEGUN 
DITORS, PUBLISHED BY THE NATIONAL BUREAU OF STANDARDS (REF. C53). THIS 
fPROXIMATION IS ACCURATE TO ATLEAST TO 7 DECIMAL PLACES. THE SYNTAX 
I3R THE FUNCTION IS: NCDF Z WHERE Z IS AN INCREASING ARRAY 

F NUMBERS FOR WHICH THE C.D.F. IS TO BE COMPUTED. THIS FUNCTION IS 
WOKED BY THE BVN FUNCTION TO COMPUTE MARGINAL PROBABILITIES. 
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V RHO BVN W;X;Y;RBAR;R;A;M;B;E;Q;D;I;J;K;Ll;L2;Al;A2;ril;M2;Pl;P2 



cn ' ' 

Z22 ‘ ' 

C33 flTHIS FUNCTION COMPUTES THE CUMULATIVE PROBABILITIES OF A STANDARD 

LAI aBIVARIATE normal DISTRIBUTION ( means 0 and st.devs 1) WITH 

C53 fiCORRELATION COEFFICIENT RHO. THE SYNTAX FOR THE FUNCTION ISi 

L61 flRHO BVN W , WHERE W=(x,y). THE FUNCTION FIRST COMPUTES A 

C73 fl CONTINGENCY TABLE OF CELL PROBABILITIES OVER A USER DEFINED PART IT 10 

LSI flOF THE X AND Y AXES AND THEN CUMULATES THE PROBAILITIES IN THE 

Z92 fi APPROPRIATE CELLS. THE USER IS PROMPTED TO INPUT THE LENGTH OF THE 

CIO] flSUBINTERVALS IN THE PARTITION OF (“4,4) (e.g., .05) AND THE DESIRED 

Cll] ft ACCURACY-^ (e.g., .0005 FOR A 3-DECIMAL ACCURACY). THE FUNCTION 

C12] flOUTPUTS BOTH PrCX<x ,Y<y3 AND PrCX>x ,Y>y3. 

C133 ' INPUT THE DESIRED PARTITION SUB INTERVAL LENGTH FOR X AND Y AXES' 

C14] 

C153 ' ' 

C16] 'INPUT THE DESIRED ACCURACY OF COMPUTATIONS' 

C17] E<-0 

CIS] X^WCl] 

C19] Y^WC2] 

C20] ■►ENDlxi. (XS“4 )vYS“ 4 
C21] ->END2x\ (X>4)'^Y>4 
C22] ^END3 xl <X>4)vy>4 
C23] X<-L/4,X 

C24] Y^L/4,Y 

C 25 ] RBAR^RHOx 1 - ( K*2 ) -5- 1 2 
C 26 ] R<-RB AR^ < 1 -RB AR*2 ) 

C27] A-«-“4+0,Kxv8-5-K 

C28] A1+(<X>A)/A) ,X, <X<A)/A 

C29] A2^( (Y>A) /A) ,Y, (Y<A)/A 

C30] Ml^( (l + Al)+“l + Al)-i-2 

C31] M2^( (l + A2)+“l + A2)-5-2 

C32] Ll^pMl 

C33] L2^/t>M2 

C34] B<-*RxMlo.xM2 

C35] P1^(NCDF l + AD-NCDF “1 + Al 

C36] P2^(NCDF 1+A2)-NCDF “14A2 

C37] REPEAT: Q^Pl-s-+/B 

C38] D<-( (LI ,L1)AQ)x( (H ,Ll)Al,Ll/>0) 

C39] B^D+.xB 
C40] Q+P2++7‘B 



C4n 

C42D 

C433 

C44: 

C453 

C46: 

C47: 

C48: 

C49D 

C503 

C513 

I 

:s2: 

:s3] 

I 

:s4: 

I 

•553 

[563 

:573 



D^< <L2,L2)|t>Q)x( (L2,L2)pl ,L2|t>0) 

B^B+.xD 

■>REPEATx^ ((+/(! P1-+/B) >E) >0) v ( + / ( | P2-+/B) >E) >0 

I<-+/X>A1 

J^+/Y>A2 

'Pr C X < ',<9X),' , Y < ',(9Y),' 3 =» 

'Pr C X > ',(5X),' , Y > ',(fY),' 3 » 

^0 

ENDl:' Pr C X < *,(»X),' , Y < ' , (?Y) , ' 

' Pr C X > ' , <5X) , ' , Y > ' , (?Y) , ' 3 

■>0 

END2: ' Pr C X < ' , (9X) , ' , Y < ' , (?Y) , ' 

' Pr C X > ',<5X),' , Y > ',(5Y),‘ 3 

^0 

END3: ' Pr C X < ',(?X),' , Y < ' , <5Y) , 

' Pr C X > ' , (?X) , ' , Y > ' , (?Y) , ' 3 



' ,5+/+/ (I , J)tB 
' ,?+/+/(I , J)+B 



1 ' 



0 ' 



1 ' 



= 0' 



= ',5<NCDF X)xnCDF Y 
',5(1-NCDF X)X(1-NCDF Y) 
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7 R'^'NCDF Z;B;p;z;T;N;W;P 

:i] flTHIS FUNCTION COMPUTES THE C. D. F. Pr C Z < z I OF A STANDARD 
C2D flNORMAL DISTRIBUTION USING THE APPROXIMATION IN EQUATION 26.2.17, 
C33 flPAGE 932 IN THE HANDBOOK OF MATHEMATICAL FUNCTIONS, EDITED BY 
ZA2 AABRAMONITZ AND STATGUN, NATIONAL BUREAU OF STANDARDS. 

C53 Z<-,Z 

Z61 l^lZAZl 
Z72 N-«-p(Z<0)/Z 

C8D M<-p<Z>0)/Z 
C9D Z^IZ 
C103 p-<-0. 2316419 

Cli: T^4-l+pxZ 

C12: z^(*-(Z*2)^-2)-5-<o2)*0.5 

C133 B-«- 0.31938153 “0.356563782 1.781477937 "1.821255978 1.330274429 

C143 P«•l-zx(To.*^5)+.xB 

C 15D ( 1-N+P) , (-M) tP 
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